library(sf)
## Linking to GEOS 3.6.2, GDAL 2.2.3, PROJ 4.9.3
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
library(leaflet)
library(scales)
library(ggmap)
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.
library(readr)
## 
## Attaching package: 'readr'
## The following object is masked from 'package:scales':
## 
##     col_factor
## Read in shapefile using sf
ak_regions <- read_sf("shapefiles/ak_regions_simp.shp")

plot(ak_regions)  

class(ak_regions)
## [1] "sf"         "tbl_df"     "tbl"        "data.frame"
head(ak_regions)
## Simple feature collection with 6 features and 3 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -179.2296 ymin: 51.15702 xmax: 179.8567 ymax: 71.43957
## CRS:            4326
## # A tibble: 6 x 4
##   region_id region      mgmt_area                                       geometry
##       <int> <chr>           <dbl>                             <MULTIPOLYGON [°]>
## 1         1 Aleutian I…         3 (((-171.1345 52.44974, -171.1686 52.41744, -1…
## 2         2 Arctic              4 (((-139.9552 68.70597, -139.9893 68.70516, -1…
## 3         3 Bristol Bay         3 (((-159.8745 58.62778, -159.8654 58.61376, -1…
## 4         4 Chignik             3 (((-155.8282 55.84638, -155.8049 55.86557, -1…
## 5         5 Copper Riv…         2 (((-143.8874 59.93931, -143.9165 59.94034, -1…
## 6         6 Kodiak              3 (((-151.9997 58.83077, -152.0358 58.82714, -1…
st_crs(ak_regions)
## Coordinate Reference System:
##   User input: 4326 
##   wkt:
## GEOGCS["GCS_WGS_1984",
##     DATUM["WGS_1984",
##         SPHEROID["WGS_84",6378137,298.257223563]],
##     PRIMEM["Greenwich",0],
##     UNIT["Degree",0.017453292519943295],
##     AUTHORITY["EPSG","4326"]]
ak_regions_3338 <- ak_regions %>%
  st_transform(crs = 3338)

st_crs(ak_regions_3338)
## Coordinate Reference System:
##   User input: EPSG:3338 
##   wkt:
## PROJCS["NAD83 / Alaska Albers",
##     GEOGCS["NAD83",
##         DATUM["North_American_Datum_1983",
##             SPHEROID["GRS 1980",6378137,298.257222101,
##                 AUTHORITY["EPSG","7019"]],
##             TOWGS84[0,0,0,0,0,0,0],
##             AUTHORITY["EPSG","6269"]],
##         PRIMEM["Greenwich",0,
##             AUTHORITY["EPSG","8901"]],
##         UNIT["degree",0.0174532925199433,
##             AUTHORITY["EPSG","9122"]],
##         AUTHORITY["EPSG","4269"]],
##     PROJECTION["Albers_Conic_Equal_Area"],
##     PARAMETER["standard_parallel_1",55],
##     PARAMETER["standard_parallel_2",65],
##     PARAMETER["latitude_of_center",50],
##     PARAMETER["longitude_of_center",-154],
##     PARAMETER["false_easting",0],
##     PARAMETER["false_northing",0],
##     UNIT["metre",1,
##         AUTHORITY["EPSG","9001"]],
##     AXIS["X",EAST],
##     AXIS["Y",NORTH],
##     AUTHORITY["EPSG","3338"]]
plot(ak_regions_3338)

ak_regions_3338 %>%
  select(region)
## Simple feature collection with 13 features and 1 field
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -2175328 ymin: 405653.9 xmax: 1579226 ymax: 2383770
## CRS:            EPSG:3338
## # A tibble: 13 x 2
##    region                                                               geometry
##    <chr>                                                      <MULTIPOLYGON [m]>
##  1 Aleutian Islands (((-1156666 420855.1, -1159837 417990.3, -1161898 416944.4,…
##  2 Arctic           (((571289.9 2143072, 569941.5 2142691, 569158.2 2142146, 56…
##  3 Bristol Bay      (((-339688.6 973904.9, -339302 972297.3, -339229.2 971037.4…
##  4 Chignik          (((-114381.9 649966.8, -112866.8 652065.8, -108836.8 654303…
##  5 Copper River     (((561012 1148301, 559393.7 1148169, 557797.7 1148492, 5559…
##  6 Kodiak           (((115112.5 983293, 113051.3 982825.9, 110801.3 983211.6, 1…
##  7 Kotzebue         (((-678815.3 1819519, -677555.2 1820698, -675557.8 1821561,…
##  8 Kuskokwim        (((-1030125 1281198, -1029858 1282333, -1028980 1284032, -1…
##  9 Cook Inlet       (((35214.98 1002457, 36660.3 1002038, 36953.11 1001186, 367…
## 10 Norton Sound     (((-848357 1636692, -846510 1635203, -840513.7 1632225, -83…
## 11 Prince William … (((426007.1 1087250, 426562.5 1088591, 427711.6 1089991, 42…
## 12 Southeast        (((1287777 744574.1, 1290183 745970.8, 1292940 746262.7, 12…
## 13 Yukon            (((-375318 1473998, -373723.9 1473487, -373064.8 1473930, -…
ak_regions_3338 %>%
  filter(region == "Southeast")
## Simple feature collection with 1 feature and 3 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 559475.7 ymin: 722450 xmax: 1579226 ymax: 1410576
## CRS:            EPSG:3338
## # A tibble: 1 x 4
##   region_id region   mgmt_area                                          geometry
## *     <int> <chr>        <dbl>                                <MULTIPOLYGON [m]>
## 1        12 Southea…         1 (((1287777 744574.1, 1290183 745970.8, 1292940 7…
pop <- read.csv("shapefiles/alaska_population.csv")
pop_4326 <- st_as_sf(pop, 
                  coords = c('lng', 'lat'),
                  crs = 4326,
                  remove = F)

head(pop_4326)
## Simple feature collection with 6 features and 5 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -176.6581 ymin: 51.88 xmax: -154.1703 ymax: 62.68889
## CRS:            EPSG:4326
##   year     city      lat       lng population                   geometry
## 1 2015     Adak 51.88000 -176.6581        122    POINT (-176.6581 51.88)
## 2 2015   Akhiok 56.94556 -154.1703         84 POINT (-154.1703 56.94556)
## 3 2015 Akiachak 60.90944 -161.4314        562 POINT (-161.4314 60.90944)
## 4 2015    Akiak 60.91222 -161.2139        399 POINT (-161.2139 60.91222)
## 5 2015   Akutan 54.13556 -165.7731        899 POINT (-165.7731 54.13556)
## 6 2015 Alakanuk 62.68889 -164.6153        777 POINT (-164.6153 62.68889)
pop_3338 <- st_transform(pop_4326, crs = 3338)
pop_joined <- st_join(pop_3338, ak_regions_3338, join = st_within)
#used st_within because we have points within polygons

head(pop_joined)
## Simple feature collection with 6 features and 8 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -1537925 ymin: 472627.8 xmax: -10340.71 ymax: 1456223
## CRS:            EPSG:3338
##   year     city      lat       lng population region_id           region
## 1 2015     Adak 51.88000 -176.6581        122         1 Aleutian Islands
## 2 2015   Akhiok 56.94556 -154.1703         84         6           Kodiak
## 3 2015 Akiachak 60.90944 -161.4314        562         8        Kuskokwim
## 4 2015    Akiak 60.91222 -161.2139        399         8        Kuskokwim
## 5 2015   Akutan 54.13556 -165.7731        899         1 Aleutian Islands
## 6 2015 Alakanuk 62.68889 -164.6153        777        13            Yukon
##   mgmt_area                   geometry
## 1         3  POINT (-1537925 472627.8)
## 2         3 POINT (-10340.71 770998.4)
## 3         4  POINT (-400885.5 1236460)
## 4         4  POINT (-389165.7 1235475)
## 5         3 POINT (-766425.7 526057.8)
## 6         4  POINT (-539724.9 1456223)
pop_region <- pop_joined %>% 
  as.data.frame() %>% 
  group_by(region) %>% 
  summarise(total_pop = sum(population), .groups = 'drop')

head(pop_region)
## # A tibble: 6 x 2
##   region           total_pop
##   <chr>                <int>
## 1 Aleutian Islands      8840
## 2 Arctic                8419
## 3 Bristol Bay           6947
## 4 Chignik                311
## 5 Cook Inlet          408254
## 6 Copper River          2294
pop_region_3338 <- left_join(ak_regions_3338, pop_region)
## Joining, by = "region"
#plot to check

plot(pop_region_3338)

plot(pop_region_3338["total_pop"])

pop_mgmt_3338 <- pop_region_3338 %>% 
  group_by(mgmt_area) %>% 
  summarize(total_pop = sum(total_pop))
## `summarise()` ungrouping output (override with `.groups` argument)
plot(pop_mgmt_3338)

plot(pop_mgmt_3338["total_pop"])

Maps

Finalizing maps with ggplot

ggplot(pop_region_3338) +
  geom_sf(aes(fill = total_pop)) +
  theme_bw() +
  labs(fill = "Total Population") +
  scale_fill_continuous(low = "khaki", high =  "firebrick", labels = comma)

# put aes in the geom_sf
rivers_3338 <- read_sf("shapefiles/ak_rivers_simp.shp")
st_crs(rivers_3338)
## Coordinate Reference System:
##   No user input
##   wkt:
## PROJCS["Albers",
##     GEOGCS["GCS_GRS 1980(IUGG, 1980)",
##         DATUM["unknown",
##             SPHEROID["GRS80",6378137,298.257222101]],
##         PRIMEM["Greenwich",0],
##         UNIT["Degree",0.017453292519943295]],
##     PROJECTION["Albers_Conic_Equal_Area"],
##     PARAMETER["standard_parallel_1",55],
##     PARAMETER["standard_parallel_2",65],
##     PARAMETER["latitude_of_center",50],
##     PARAMETER["longitude_of_center",-154],
##     PARAMETER["false_easting",0],
##     PARAMETER["false_northing",0],
##     UNIT["Meter",1]]
ggplot() +
  geom_sf(data = pop_region_3338, aes(fill = total_pop)) +
  geom_sf(data = rivers_3338, aes(size = StrOrder), color = "black") +
  geom_sf(data = pop_3338, aes(), size = .5) +
  scale_size(range = c(0.01, 0.2), guide = F) +
  theme_bw() +
  labs(fill = "Total Population") +
  scale_fill_continuous(low = "khaki", high =  "firebrick", labels = comma) +
  scale_x_continuous(breaks = seq(-180, 180, by = 20))

Raster basemaps

pop_3857 <- pop_3338 %>%
  st_transform(crs = 3857)
# Define a function to fix the bbox to be in EPSG:3857
# See https://github.com/dkahle/ggmap/issues/160#issuecomment-397055208
ggmap_bbox_to_3857 <- function(map) {
  if (!inherits(map, "ggmap")) stop("map must be a ggmap object")
  # Extract the bounding box (in lat/lon) from the ggmap to a numeric vector, 
  # and set the names to what sf::st_bbox expects:
  map_bbox <- setNames(unlist(attr(map, "bb")), 
                       c("ymin", "xmin", "ymax", "xmax"))
  
  # Coonvert the bbox to an sf polygon, transform it to 3857, 
  # and convert back to a bbox (convoluted, but it works)
  bbox_3857 <- st_bbox(st_transform(st_as_sfc(st_bbox(map_bbox, crs = 4326)), 3857))
  
  # Overwrite the bbox of the ggmap object with the transformed coordinates 
  attr(map, "bb")$ll.lat <- bbox_3857["ymin"]
  attr(map, "bb")$ll.lon <- bbox_3857["xmin"]
  attr(map, "bb")$ur.lat <- bbox_3857["ymax"]
  attr(map, "bb")$ur.lon <- bbox_3857["xmax"]
  map
}
bbox <- c(-170, 52, -130, 64)   # This is roughly southern Alaska
ak_map <- get_stamenmap(bbox, zoom = 4)
## Source : http://tile.stamen.com/terrain/4/0/4.png
## Source : http://tile.stamen.com/terrain/4/1/4.png
## Source : http://tile.stamen.com/terrain/4/2/4.png
## Source : http://tile.stamen.com/terrain/4/0/5.png
## Source : http://tile.stamen.com/terrain/4/1/5.png
## Source : http://tile.stamen.com/terrain/4/2/5.png
ak_map_3857 <- ggmap_bbox_to_3857(ak_map)

ggmap(ak_map_3857) + 
  geom_sf(data = pop_3857, aes(color = population), inherit.aes = F) +
  scale_color_continuous(low = "khaki", high =  "firebrick", labels = comma)
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

  # scale_x_continuous(breaks = seq(-170, 130, by = 10))

Leaflet

epsg3338 <- leaflet::leafletCRS(
  crsClass = "L.Proj.CRS",
  code = "EPSG:3338",
  proj4def =  "+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs",
  resolutions = 2^(16:7))

st_crs(pop_region_3338)
## Coordinate Reference System:
##   User input: EPSG:3338 
##   wkt:
## PROJCS["NAD83 / Alaska Albers",
##     GEOGCS["NAD83",
##         DATUM["North_American_Datum_1983",
##             SPHEROID["GRS 1980",6378137,298.257222101,
##                 AUTHORITY["EPSG","7019"]],
##             TOWGS84[0,0,0,0,0,0,0],
##             AUTHORITY["EPSG","6269"]],
##         PRIMEM["Greenwich",0,
##             AUTHORITY["EPSG","8901"]],
##         UNIT["degree",0.0174532925199433,
##             AUTHORITY["EPSG","9122"]],
##         AUTHORITY["EPSG","4269"]],
##     PROJECTION["Albers_Conic_Equal_Area"],
##     PARAMETER["standard_parallel_1",55],
##     PARAMETER["standard_parallel_2",65],
##     PARAMETER["latitude_of_center",50],
##     PARAMETER["longitude_of_center",-154],
##     PARAMETER["false_easting",0],
##     PARAMETER["false_northing",0],
##     UNIT["metre",1,
##         AUTHORITY["EPSG","9001"]],
##     AXIS["X",EAST],
##     AXIS["Y",NORTH],
##     AUTHORITY["EPSG","3338"]]
pop_region_4326 <- pop_region_3338 %>% st_transform(crs = 4326)

pal <- colorNumeric(palette = "Reds", domain = pop_region_4326$total_pop)

m <- leaflet(options = leafletOptions(crs = epsg3338)) %>%
        addPolygons(data = pop_region_4326, 
                    fillColor = ~pal(total_pop),
                    weight = 1,
                    color = "black",
                    fillOpacity = 1,
                    label = ~region) %>% 
        addLegend(position = "bottomleft",
                  pal = pal,
                  values = range(pop_region_4326$total_pop),
                  title = "Total Population")

m